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prior to its application for cell design or optimisation purposes. In this fitting, a set of experimental data 
is used to guess the value of those parameters of the model that cannot be either modelled or measured 
experimentally. This is crucially the case of the charge transfer coefficients (ap, đa, @b,c» @tc) and the 
exchange current densities (ip, io.) in the Butler-Volmer equation (i.e. electrochemical model). 

Keywords: The fitting of the electrochemical parameters in the SOFC, SOEC and SORFC modelling literature is 
SOFC reviewed in this work. It is found that this process is only vaguely discussed, if mentioned at all. In the 
a er fitting authors’ opinion, this practice contributes with uncertainty rather than guidance, since this fitting process is 
Exchange current density of utmost significance for making reliable quantitative predictions. i , , 
Transfer coefficient In this work, we further introduce a comprehensive model for the simulation of solid oxide regenerative 
Modelling fuel cells, ie. a model that simulates, without any ad hoc adjustments or tuning, both the SOFC and SOEC 
modes. We also describe in detail how the electrochemical parameters are fitted, and discuss the applicability 
of the values commonly used in the literature for these fitted parameters and their proper validation. Finally, 
the validity of the proposed model and fitted parameters is shown by comparison of the numerical results 
with experimental data. 
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1. Introduction 


Hydrogen has been identified as a promising energy carrier: it 
can be cleanly produced in a water electrolyser driven by renew- 
able energy; it may be then stored or transported, and eventually 
converted back into water and electricity in a fuel cell. A Solid 
Oxide Regenerative Fuel Cell (SORFC) is an electrochemical device 
that allows both operating modes, i.e. it can operate reversibly 
either as a fuel cell or as an electrolyser. 

An SORFC operating as a solid oxide fuel cell (SOFC) generates 
electricity, heat and water from hydrogen and oxygen. An SORFC 
operating as a solid oxide electrolyser cell (SOEC) produces hydrogen 
and oxygen when water, heat and electricity are supplied. 

The reversibility of the solid oxide cells was first proved in the 
1990s for both tubular and planar geometries [3-5]; however, 
research in this field subsided due to the low prices of fossil fuels. 
In recent years environmental, economic and geopolitical concerns 
over fossil fuels have rekindled the interest in this technology. This 
is evinced by the ongoing research work for both planar [6-8] and 
microtubular solid oxide regenerative fuel cells [9,10]. 

Numerical modelling is a primary tool to understand, and 
optimise, the operation of solid oxide cells, either in fuel-cell or 
electrolyser mode. The modelling of either SOFCs or SOECs has 
been the subject of many papers [11-19]. Often such models use 
computational fluid dynamics (CFD) to simulate all the relevant 
mass, heat and charge transfer processes. Wang et al. [1] present 
an overview of the several modelling alternatives, including: 
physical models (by which they mean those that represent, 
mathematically, the underlying physics, be it in zero, one, multiple 
spatial dimensions); equivalent circuit models (based on electro- 
chemical impedance spectroscopy, or EIS, measurements); and 
“gray-” and “black-box” models (such as Artificial Neural Networks 


Table 1 


and neuro-fuzzy systems). The review paper by Hajimolana et al. 
[2] is a systematic inventory of the main submodels proposed 
for the mathematical representation of all the relevant physical 
processes in the several spatial domains (gas channels, electrodes, 
electrolyte, interconnects). 

References reporting the modelling of both SOFC and SOEC 
operating modes (i.e. SORFC) with a single CFD model that works 
satisfactorily in both situations are very scarce. To the authors’ 
knowledge only Jin and Xue [20] have presented a validated 
numerical tool for the simulation of SORFCs. (Ni et al. [21,22] 
address the modelling and validation of a reversible solid oxide 
fuel cell, but in fact they only report the SOEC behaviour.) If 
properly formulated, a CFD model for either SOFCs or SOECs 
should produce suitable results in the other regime without any 
modification to the laws (submodels) for the underlying funda- 
mental physics. However, this is not the case for most models, as 
shown below. This paper reviews the challenges posed by the 
development of such a unified model, and how to overcome them. 

Any comprehensive SOFC, SOEC or SORFC model relies to a 
certain extent on data fitting to find values for some of the physical 
parameters, as it is the case of the charge transfer coefficients 
(dba, Cfa,%b,c, Mc) and the exchange current densities (io.a, io) in 
the electrochemical model. This fitting, when properly resorted to, 
is the consequence of the incomplete knowledge of the underlying 
physical phenomena, or of the excessive complexity of such 
phenomena for them to be accommodated within a fluid-flow 
model. The disparity of spatial scales between these phenomena 
and the device to be simulated is often one of the sources of this 
complexity. 

However, often this fitting process is only vaguely discussed 
in the literature. Table 1 summarises the common practices in 
the SOFC, SOEC and SORFC modelling literature to evaluate the 


Review of the evaluation of the electrochemical parameters in SOFC, SOEC and SORFC modelling literature,f(©) means “calculated as a function of ©”. 


Model Fitting Opa / Aba ap,c/ Ac 
SOFC 

[25] No 0.5/0.5 0.5/0.5 
[57] No 0.5/0.5 0.5/0.5 
[50] No 0.5/0.5 0.5/0.5 
[12] No 0.5/0.5 0.5/0.5 
[59] No 0.5/0.5 0.5/0.5 
[60,15] Yes 0.5/0.5 0.5/0.5 
[16] Yes 0.5/0.5 0.5/0.5 
[45] Yes 1.5/0.5 0.75/0.25 
[28] Yes 2/1 [23] 1.5/0.5 [24] 
SOEC 

[61,62] No 0.5/0.5 0.5/0.5 
[29] Yes 0.5/0.5 0.5/0.5 
[63] No - - 
[21,17,26] No 0.5/0.5 0.5/0.5 
[22] Yes 0.5/0.5 0.5/0.5 
[19] Yes 0.5/0.5 0.5/0.5 
[51] No 0.5/0.5 0.5/0.5 
SORFC 

[20] Yes 0.5/0.5 0.5/0.5 


loa [ioe 


f(T) 
5300/2000 A m~? 
5300/2300 A m~? 


7460/10090 A a 


5300/2000 A m~? 
Both f(T, species) 
Both f(T, species) 
Both f(T, species) 
Both f(T, species) 
Both from [25] 
Kinetics/fitted 
2000/5300 A m~? 
200 Am~?/ 

Fitted 

Both f(T, species) 


1/01 A m~? 


Value source 


From [56]/ guessed 


From [58] 
Not justified 
Not justified 
From [57] 
Fitting 
Fitting 
Fitting 
Fitting 


SOFC 

Poor fitting 

Not mentioned 
From [57] 

Not mentioned 
Poor agreement 
Literature 


Fitting 
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Nomenclature 

a exponent, Eq. (52), dimensionless 

A constant in the ionic conductivity equation, Eq. (A.19), 
AV-'m-! 

ie JANAF constant for species a, dimensionless 

Ary JANAF constant for species a, K7! 

aza JANAF constant for species a, K7? 

Aa JANAF constant for species a, K~? 

sq JANAF constant for species a, K~* 

b exponent, Eq. (52), dimensionless 

B constant in the ionic conductivity equation, 
Eq. (A.19), K 

Bo porous-medium permeability-coefficient, m? 

Ç exponent, Eq. (53), dimensionless 

Ca molar density of species a, kmol m~? 

Cp specific heat at constant pressure of the fluid, 
m? s~? K7! 

Cpa specific heat at constant pressure of species a, 
m°? s7? KT! 

dp mean pore diameter, m 

Dap ordinary diffusion coefficient of species a in species 2, 
ms! 

Dam ordinary diffusion coefficient of species a in a 
multicomponent mixture, m? s7! 

Dru Knudsen-diffusion coefficient of species a, m? s7! 

Enact activation energy, kg m? s7? kmol~! 

E, black-body emissive power, kg s~? 

E electric potential difference, V 

EY standard electric-potential difference, V 

F Faraday's constant, As kmol~! 

Fi_j view factor between the i-th and the j-th surface 
elements, dimensionless 

h sensible enthalpy of the fluid, m? s7? 

ha sensible enthalpy of species a, m? s~? 

Ho external irradiation, kg s~? 

i, current density magnitude, A m~? 

i current density, A m7? 

io, exchange current density, A m7? 

Ja mass diffusion-flux of species a relative to the 

a mass-average velocity, kg m? s~! 

ome, superficial mass diffusion-flux of species a through a 
porous media, kg m° s~! 

kg Boltzmann constant, kg m? s7? K7! 

n number of electrons released during the ionisation 

= process of one fuel molecule, dimensionless 

Na total molar flux of species a, kmol m~? s~! 

p pressure, kg m~! s7? 

Pa partial pressure of species a, kg m7! s7? 

Ci energy flux, kg s~? 

Q volumetric heat sources/sinks, Joule heat, kg m~! s7? 

r length from the i-th to the j-th surface elements, m 

R ideal gas constant, kg m? s72 kmol ~! K~! 

Sy, mass source or sink of species a, kg m? s~! 

Soia Sutherland-law parameter for species a, K 

Sona Sutherland-law parameter for species a, K 

Sx, molar source or sink of species a, kmol m? s~! 

t time, s 

T temperature, K 

T parameter defined in Eq. (A.7), dimensionless 

Toia Sutherland-law parameter for species a, K 

Toya Sutherland-law parameter for species a, K 

v fluid mass-averaged velocity, m s~! 


superficial permeation velocity, m s~! 
voltage, V 

molecular weight of the fluid, kg kmol~! 
molecular weight of species a, kg kmol”! 
molar fraction of species a, dimensionless 
mass fraction of species a, dimensionless 


Greek letters 


a species, dimensionless 

db backward transfer coefficient, dimensionless 

a forward transfer coefficient, dimensionless 

y pre-exponential coefficient, A m~? 

Fy dusty-gas model parameter, Eq. (12), m? s~! 

i dusty-gas model parameter, Eq. (12), kg~! s! kmol 

AH; specific enthalpy of reaction, kg m? s72? kmol7! 

e porosity of the porous medium, dimensionless 

Ëx characteristic energy of species a , kg m? s~? 

Nact activation overpotential, V 

Ncon concentration overpotential, V 

fohi Ohmic overpotential, V 

0 angle, 

A thermal conductivity of the fluid, kg m s73 K7! 

An thermal conductivity of species a, kg m s7? K7! 

Joa Sutherland-law parameter for species a, kg m s~? K7! 

u fluid viscosity, kg m~! s~! 

Ha viscosity of species a, kg m~! s7! 

ie Sutherland-law parameter for species a, kg m~! s~! 

v number of occurrences of the rate determining step in 
an overall reaction, dimensionless 

é surface emissivity, dimensionless 

p mass density of the fluid, kg m~? 

Pp charge density, C m~? 

č ionic conductivity, A V7! m~! 

Oa characteristic length of species a, A 

Cap average collision diameter of species a and £, A 

Osh Stefan-Boltzmann constant, kg s73 K~4 

T tortuosity factor of the porous medium, dimensionless 

vo . dusty-gas model parameter, Eq. (13), ms~! 

3 dusty-gas model parameter, Eq. (13), 
kg~' m-1 s! kmol 

ca dusty-gas model parameter, Eq. (14), ms~! 

T, dusty-gas model parameter, Eq. (14), 
kg~' m-1 s! kmol 

h electric potential, V 

pa electric potential of the electric-conducting phase of 
the fuel electrode, V 

bc electric potential of the electric-conducting phase of 
the air electrode, V 

de electric potential of the ion-conducting phase of the 
electrolyte, V 

pe effective electric potential of the ion-conducting phase 
of the electrolyte, V 

Bap semi-empirical correlation parameter, 
Eq. (A.12), dimensionless 

o heat of reaction, kg m~! s~? 

Qp collision integral, dimensionless 

Subscripts 

a fuel electrode 

b bulk, feeding mixture concentration 
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Ç air electrode 

e electrolyte 

eq equilibrium 

f fluid 

i isothermal-surface index 

ia interconnect/fuel—-electrode interface 

ic interconnect/air—electrode interface 

irrev irreversible 

J isothermal-surface index 

oc open circuit 

r generic reaction wall, i.e. electrode—-electrolyte 
interface 

ra anode reaction wall, i.e. anode-electrolyte interface 


rc cathode reaction wall, i.e. cathode-electrolyte 
interface 

rad radiation 

ref reference value 

rev reversible 

s solid 

SOEC at SOEC mode 

SOFC at SOFC mode 

Superscripts 

eff effective value in the porous media 


electrochemical parameters. It is widely assumed that the charge 
transfer coefficients are equal to 0.5, what would imply a single step, 
single electron reaction. Theoretically, it is clear that this assumption 
does not apply to the hydrogen reduction/oxidation reaction. Experi- 
mentally, there are interesting works refuting this hypothesis [23,24]. 
Numerically, there is also an interesting contribution of Grondin et al. 
[19] stating that Butler-Volmer equation together with the single 
step, single electron hypothesis does not permit relevant computing 
predictions of SOEC polarisation curves. However, this assumption is 
still used as it is simple and provides (apparently) good fitting to 
experimental data. In this work, we show that this assumption fails if 
the model is properly validated and that more accurate values can be 
found by means of a thorough, careful fitting process. 

Once the electrochemical reaction is assumed as a single step, 
single electron one, the exchange current densities are the 
remaining unknowns of the model. How these are found and 
reported in the published literature deserves some close scrutiny. 
As shown in Table 1, some papers in the literature do not mention 
these values. Sometimes, when no experimental data is found 
these values are guessed without any further validation [25]. Some 
authors take constant values from published experimental works, 
even if the operating conditions or materials of the experimental 
cell do not match those of the cell being modelled [21,26]. For 
instance, Ni et al. [26] assume that the exchange current density of 
the carbon dioxide electrolysis is that of the water electrolysis. Ni 
et al. [21] validate the numerical data with the experimental data 
of Momma et al. [27], at first sight this validation would indicate 
that taking constant values from the literature may lead to 
reasonable fitting with experimental data. However, close exam- 
ination of the experimental data used [27] reveals that only the 
linear range of the experimental characteristic curve is used 
(i e (0,5000) Am~?), avoiding the comparison with the experi- 
mental data at high current densities where the experimental 
trend is not linear anymore. Table 1 also shows that the value of 
the exchange current densities is estimated, in other works, to fit 
the numerical results to an experimental data set. For example, in 
[20] a constant value for both exchange current densities is 
guessed to fit the experimental data, neglecting its dependence 
on temperature and species concentration; the resulting values are 
unusually small (1,0.1 Am~?). This is also done in [22] for one 
exchange current density, but the authors do not indicate whether 
it is the cathodic or the anodic one. In the SOFC literature, it is a 
common practice to use experimental correlations that express the 
exchange current densities as a function of the temperature and 
species distribution at the electrodes, e.g. [28]. These are generally 
Arrhenius type expressions, the pre-exponential factors of which 
are the fitted parameters. In the SOEC literature this kind of 
correlation is used in [29] leading to a very poor agreement. 

From the above review, the authors conclude that the evalua- 
tion of the electrochemical parameters is, in general, not properly 


justified. In our opinion, this not only greatly diminishes the 
usefulness of modelling as a design tool, but also contributes with 
uncertainty rather than guidance when published in the literature. 

In the present work we introduce a comprehensive model for 
the simulation of an SORFC, together with a detailed description of 
how the electrochemical parameters are fitted. The validity of the 
proposed model and fitted parameters is also shown by compar- 
ison of the numerical results with experimental data for a micro- 
tubular SORFC [10]. 

The model presented caters for the following phenomena: 
(a) mass transport through channels; (b) momentum transport 
through channels and porous solids; (c) multispecies mass transfer 
through porous media (convection, ordinary diffusion, Knudsen 
diffusion); (d) heat transfer in gases, porous media and solids 
by conduction, convection and radiation; (e) charge transport 
through an impervious solid; and (f) reversible electrochemistry. 
The model has been implemented in OpenFOAM-1.5-dev [30], an 
open-source CFD platform based on the finite-volume method. 
The full mathematical model and the relevant features of the 
numerical method are thoroughly described in Sections 2 and 3 
respectively. The model can be regarded as a compendium of the 
most appropriate submodels found in the literature for the task; 
and the full description of this comprehensive model in this paper 
should enable its implementation in any code by any interested 
researcher. 

The model is validated using experimental data presented in 
Section 4. The unavoidable fitting process is next described in 
Section 5. It reveals that some common simplifications used in 
the SOFC and SOEC modelling literature (see Table 1) might be 
inappropriate. Our findings highlight the often-ignored relevance 
of this fitting process for obtaining a predictive tool that can be 
used reliably across the range of SOFC, SOEC and SORFC designs. 

Finally in Section 6 our numerical results are compared with 
experimental data to validate the model. Not only the numerical 
and experimental characteristic i-V curves agree well, but so do 
the magnitudes of other variables that depend on the fitted 
parameters, such as the exchange current density and the activa- 
tion overpotential. 


2. Mathematical model 


In this section the mathematical model that describes the 
operating principles of an SORFC is presented. The model is an 
extension of one developed previously for SOFC simulation [16,31], 
the new features being: (i) a new multidimensional charge- 
transfer model for the electrolyte; (ii) an improved formulation 
of the calculation of heat sources and sinks; (iii) the extension of 
the electrochemical model to SOEC conditions. 
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The operating principle of a solid oxide regenerative fuel cell, both 
in SOFC and SOEC modes, involves a complex network of interde- 
pendent physical phenomena within the several SORFC components. 
The governing equations that describe a given physical phenomenon 
may change from one component to another because of the 
components’ differing nature, e.g. porous media, gas channels or 
impervious solids. To describe each phenomenon with the appro- 
priate equations, the mathematical model has been split into a set of 
submodels, namely: (i) channel; (ii) electrode; and (iii) electrolyte 
submodels. Two additional submodels account for the radiative heat 
transfer and the electrochemistry; these operate at the interface 
between some of the components (respectively channel-electrode 
and electrode-electrolyte). 


2.1. Channel model 


The SORFC channels are ducts along which the reactant 
mixtures are driven from the SORFC inlet to the electrodes. The 
byproducts of the electrochemical reactions exit the electrode and 
join the unconverted reactants in the channel, where they are 
evacuated to the SORFC outlet. The channel model accounts for the 
phenomena summarised in Table 2 and described below. 

The differential form of the continuity equation is 


op —. = 
StV-(p¥)=0 a) 


where V is the fluid velocity and p is the fluid density, given by 
the ideal gas law: 


Ww 


P= RT 2) 


where p is the pressure, W the mixture molecular weight, T the 
temperature and R the ideal gas constant. 

The momentum-conservation equation for a Newtonian fluid 
with negligible bulk viscosity (wy + 0), with negligible body forces, 
is (in vector form) 


Apv 
we Jiv.@vv¥ 


==> 
J-V-(VV)=V- MITY -uvv T -vp 


(3) 


where y is the dynamic viscosity of the fluid (see Appendix A.1). 
The differential form of the equation for the conservation of 

chemical species a in the channel is 

ol > 

Me) Y. oy, V): Ta =S), (4) 

where y, is the mass fraction of species a, F, is the mass 


diffusion-flux of species a relative to the mass-average velocity, 
and Sy stands for the volumetric sources or sinks. 


Table 2 
Summary of the main features of the mathematical submodels. 


Submodel Phenomena State/constitutive equations 
Channel Mass transport Ideal gas 

Momentum transport Newtonian fluid 

Species transport Consistent diffusion 

Energy transport Fourier multicomponent 
Electrode Mass transport Ideal gas 

Momentum transport Darcy 

Species transport Dusty-gas 

Energy transport Fourier multicomponent 
Electrolyte Charge transport Ohm 

Energy transport Fourier 
Radiation Surface-to-surface View-factor method 


Electrochemical Redox Nernst, Butler-Volmer 


We model F, according to the consistent effective binary 
diffusion method [32]: 


-> 
J a = —PDamVYa +PYa 2 (DpmVYp) (5) 
vB 


where Dam is the diffusion coefficient of species a in a multi- 
component gas mixture (see Appendix A.2). This formulation 
ensures that the diffusive fluxes relative to averaged fluid velocity 
properly add up to zero. It thus combines the simplicity of the 
effective binary diffusion approach with flux consistency. 

From Eqs. (4) and (5), the conservation equation of a chemical 
species a in the SORFC channel is 


ol 

NYa) Y. (0y, V)—V (Dam VY) +Y [rgo vyp] =s, (6) 
VB 

The conservation of energy in the system is accounted for 


through the following equation for the sensible enthalpy: 


aT a 
PC +Y -(pCp VT)—V - AVT) = — Daher + TV 
Va 


(Caan ya [ZT aha) (7) 


where C, is the fluid specific heat at constant pressure 
(see Appendix A.5), 4 is the thermal conductivity of the fluid 
(see Appendix A.4), and h, is the sensible enthalpy of species a 
(see Appendix A.6). 

In Eq. (7) the following assumptions and provisions have been 
made [16]: (i) the sensible enthalpy is defined as dh=C,dT; 
(ii) the specific heat at constant pressure depends on temperature 
Cp =f(T) (see Appendix A.5); (iii) the diffusion of heat is due to 
both heat conduction (Fourier's law) and species mass diffusion; 
(iv) the power delivered to the fluid by body forces is negligible; 
(v) viscous dissipation of energy can be neglected; (vi) large 
pressure differences in the system are not expected; (vii) the 
reaction heat-release within the channel is null since the reaction 
sites are located in the electrodes; and (viii) radiation is from 
surface to surface, and hence there are no volumetric sources of 
radiative heat in the channels [33,34]. 


2.2. The electrode model 


The electrode model accounts for the momentum, mass, species 
and heat transport within the SORFC electrodes; both anode and 
cathode are porous ceramics that are assumed to be homogeneous 
and isotropic. 

The conservation equation for the species @ in a porous 
medium, i.e. the counterpart in the electrode to Eq. (4) in the 
channel, is 


(EPY a) 
ot 
where (V) is the fluid permeation velocity through the electrode 
(or superficial velocity) and (j ,)=eJj , represents the superficial 
combined diffusion flux (ordinary and Knudsen) through the 

porous medium. 

The dusty gas model (DGM) is the most convenient approach to 
model such a species transfer, as reported in [35,36]. It is worth 
noting, though, that the theoretical basis of the DGM [37] was 
criticised by Kerkhof [38,39], who as an alternative proposed the 
binary-friction model (BFM). However, the BFM does not improve 
the DGM results [38] and both models showed similar agreements 
and discrepancies with experimental data. The DGM is thus used 
in this work; its general form is as follows [37]: 


+V- (Cy AV)L+V Ff a) = Sya (8) 


1 Vp 5 XN a—XaN 5 Na 
RT S” pHa oe 


1 Pa Bo 
off off 


Vp (9) 
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=i 
where N, is the total molar flux of species a, Bọ is a constant 


D ines permeability-coefficient, and oy =(e/t)Dyg and 


Di. f (e /t)Dx stand for the effective binary and Knudsen diffu- 
sion coefficients respectively (see Appendices A.2 and A.3). 


Considering _the total molar flux of species a, Ng= 
W2'(py,(V¥)+(F a), Eq. (8) is rewritten as 
OlECa) 
aot Na = Sra (10) 


where Sxa is the volumetric molar sources or sinks of species a, 
Sya = W, 1Sa. 
Replacing Eq. (9) into Eq. (10) and rearranging it results in 


HE Ca) 
ot 


-V ($VP) +V: (PP PA +Y ON py) = Sr (11) 


where the following coefficients, here referred to as DGM para- 
meters, have been defined for compactness and ease of interpreta- 
tion: 


» Ig 1 


 —— = 
1 
Pa r) + ay] 


«= RT (12) 


—>P 

=p 0, sp To Bo 

Va — RT Fal vp) (13) 
—>N vid 

sN v N Np 

Ua =p > 0, =TaGRT (14) 
RT yen (Fr) 


I, (m? s~') represents a global diffusion coefficient of species a in 
the multicomponent gas mixture in a porous media ; cua (m s~!) 
is a superficial velocity due to the pressure gradient in the porous 
medium; and a. (m s7!) is a superficial velocity induced by the 
flow of the other species. 

Eq. (11), applied to all the components of the gas mixture, 
represents the complete set of equations for mass transport, and 
do not need to be supplemented by any additional equations of 
motion [37]. The global pressure p can be computed as 


p= Pa (15) 


and the superficial permeation velocity may then be estimated 
using Darcy's Law [40]: 


ty = -Pvp (16) 
H 


The authors have published an open-source species mass transfer 
library that includes this DGM algorithm, it is available in [41]. 
The sensible enthalpy conservation equation in a porous 
medium, assuming local thermal equilibrium between the porous 
matrix and the fluid flowing through the interstitial space [42], is 


Aepphy +A —e)pshs] 
at 


EV- hY- leda- (17) 


where the subscripts f and s stand for fluid and solid respectively, 
and the following assumptions and provisions have been made: 
the power delivered to the fluid by the body forces and the viscous 
dissipation are negligible; there are no large pressure differences 
in the system; and there are no volumetric heat sources within 
the electrodes because Joule heating takes place only within the 
electrolyte, and radiation can be neglected in the opaque electro- 
des [34,33]. Note that surface-to-surface radiative heat flux, jad, 
is considered at the electrode outer boundary by means of a 
boundary condition (see Section 2.4). 


The diffusive heat fluxes appearing in Eq. (17) are given by 


> Ta 
df= —4 VT+ 2 alla) (18) 
Gs=—AsVT (19) 


An effective thermal conductivity of the porous medium can be 
defined as 4f = elf +(1-— e)s [43]. Considering Eqs. (18), (19) and 
the following relationship: 


> — 
prh V + XG aha) = Liha Wa a) (20) 
then Eq. (17) can be rewritten, after some rearrangement, as 
oT 
le Cop + -9Co k+: [TCW No] -v -0V 
Va 


=TV. [E CWN o] — dlhaV . (WaN a)] 


Va 


eh i (1—eh, Ps + (21) 


where o, the heat of reaction, has two contributions: (a) due to the 
reversible reaction heat release (Qey), this being released only in 
the fuel electrode; and (b) due to the irreversibilities of the process 
(Qirrev), this being released where such irreversibilities or over- 
potentials take place. Thus, 


Qirrev,c 
Q= 
Qrev P Qirrev.a 


The chemical reaction is assumed to take place on the electrode- 
electrolyte interfaces, and thus this is where the heat of reaction is 
released; it is therefore accounted for in Eq. (21) as a heat flux 
boundary condition (see Section 3). 


in the air electrode 


in the fuel electrode (22) 


2.3. The electrolyte model 


The electrolyte is an impervious solid traversed by an anionic (O77 ) 
flux. The model consists of equations for the conservation of charge, 
which simultaneously enforces mass conservation, and of energy. 

The charge-conservation equation is, considering Ohm's Law, as 
follows: 


op 
ot 
where ò is the volumetric charge density, e is the ionic con- 
ductivity of the electrolyte (see Appendix A.7), and e is the 
electric potential of the ion-conducting phase of the electrolyte. 
Heat transfer in the electrolyte is due to conduction. Considering: 
(i) Fourier's law, d = —AVT; (ii) Joule heating, Q=i*é.~'; and 
(iii) the definition of sensible enthalpy, dh = Cp dT, then the energy 
conservation equation can be written in terms of temperature as 
op i? 
— hogt A (24) 


—V : (Ce Vpe)=0 (23) 


o- V- (AaVT)= 


where i = Ti is the current density magnitude. 

For those cell configurations where the electrolyte is exposed 
to the channels, such as the case of a microtubular cell where 
the cathode is shorter than the electrolyte, the surface-to-surface 
radiative heat exchange is included as a heat-flux boundary condition 
using a surface-to-surface radiation-model (see Section 2.4). 


2.4. The radiation model 


The importance of radiative heat transfer has been discussed in 
the SOFC-modelling literature [33,34,14]. The prevailing ideas are: 
(i) the most common gases in SOFCs behave as non-participating 
media, and thus surface-to-surface radiation exchange is the 
only radiative transfer mode that needs to be considered in the 
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channels of an SOFC; (ii) the electrodes are opaque to radiation; 
(iii) the effect of volumetric radiation within a thin electrolyte is 
negligible, though the electrolyte material is optically thin. These 
conclusions are also applicable to SOEC. 

The accurate calculation of surface-to-surface radiation requires of 
complex models accounting for both parallel and oblique radiation 
[14], such the view-factor method used in the present work [44]: 


1 
radi = aE — LFi-jEbj)—Hoi+ È i(=- 1) Fira l (25) 
j J LNS 


where qaq; is the surface-to-surface radiative heat flux arriving at the 
infinitesimal, isothermal surface element dA;; and j indicates an 
infinitesimal, isothermal element of surface dA; irradiating element i. 
In Eq. (25), £ is the emissivity of the surface, E,=0o,,T* is the 
blackbody emissive power, Gsp being the Stefan—Boltzmann constant; 
H, is the incident radiation entering or leaving the domain through an 
opening (H,=0, in this work), and F;_; is the view factor between the 
two infinitesimal elements dA; and dAj. The view factors are defined as 


cos 0; COS 0j 


Fi-j= ar? 


dAj (26) 
where 6; and 6; are the angles between the normal vector to the 
surface elements dA; and dA; and the line (of length r) connecting 
them. Further details of the radiation model are given elsewhere | 16]. 


2.5. The electrochemical model 


The electrochemical model represents the electrochemical 
reaction that takes place at the triple phase boundaries (TPBs), 
where electrolyte, reactants and the catalyst region of the elec- 
trode meet. These TPBs are confined to a region close to the 
electrode-electrolyte interface and therefore the reaction is repre- 
sented as occurring on a surface. 

In the following subsections the electrochemical models for the 
SOFC and SOEC operating modes are presented. 


2.5.1. SOFC mode 

The electrochemical reaction that takes place in a hydrogen-fed 
SOFC is 
H2 +4 O2 > H20 


The standard expression for the cell voltage of such a reacting 
system is 


V=Voc "Nconc,a — Nconc,c — ohm — Nact,a — Mact,c (27) 


where subscripts a and c stand for the fuel and air electrodes 
respectively (see Table 3) and yeonce Nonm, Nact are the concentration, 
ohmic and activation overpotentials respectively. These overpo- 
tentials, which are positive, are defined below: 


è Concentration overpotential: 


RT Pu, bPH,0,ra 
1 = In 2l 20, (28) 
ionet OR (Baeten) 


Neonc,c = = In zt (29) 
i 2F V Po,,rc 


Table 3 
Electrochemical roles of the cell electrodes for both operating modes. 


Electrode Label SOFC mode SOE mode 
Fuel electrode a Oxidation Reduction 
Air electrode € Reduction Oxidation 


© Ohmic overpotential: From Eq. (27) 


ohm = Voc — V — Nconc,a —Mconc,e — Nact.a — Nact.c (30) 
where by definition [45] 

acta = Ea — Eeq,a (31) 
Nact,c = Eeg, — Ec (32) 
Voc = Vrev+Nconca + Mconc.c (33) 
with Eega, Eeq,c, Ea, Ec and Vrey being defined as 

Pog | E? a In ( Puss) (34) 

RT- ee ae 

Eegc = JF In ( Prep ) | (35) 
Ea = ha — bea (36) 
Ec = pe — ec (37) 
Vrey = Eeq,c — Eeq,a (38) 


Here pa and c are the electric potentials of the electron- 
conducting phase in the fuel and air electrodes; and ea and 
pec are the electric potentials of the ion-conducting phase in 
the fuel and air electrodes respectively. 

From Eqs. (31) to (38), Eq. (30) may be rewritten as 


Nohm = Pera = Perc H Perc ~ Para -V (39) 


In the present work, the electrodes are assumed to be ideal 
electron conductors: 


Vøa=0 = ¢,=constant => Para = Paia (40) 


Vé-=0 =¢,=constant => Pere = Pic (41) 


where aia and ¢,j- are the electric potential of the electron- 
conductor phase of the electrode at the electrode-interconnect 
interface, that satisfy the following equation: 


V= Pcic = Paia (42) 


Therefore, from Eqs. (40), (41), (42) and (39), the ohmic over- 
voltage is defined as 


Nohm = Pera — Pere (43) 


The ohmic overpotential is, in Eq. (27), defined as a positive 
magnitude; therefore from Eq. (43): 


Nohm > 0 => Pera > Perc (44) 
Eq. (44) conflicts with Ohm's law (Qera < erc) as the current 


flows from the cathode to the anode. To remedy this, a new 
variable is used: 


be = — be (45) 

Eqs. (31) and (32) may be rewritten as 

Nact,a = Para + Pera — Eeq,a (46) 

Nacte = Eege — bere — Perc (47) 

From Eqs. (30), (46) and (47) 

Nohm = Pere — Pera + Pere — Para — V (48) 
Vv 


Therefore, the ohmic overpotential is redefined as 


Nohm = Pire a Pera (49) 
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where onm is positive as required and ¢* satisfies Ohm's law 
(Pera < Perc): 

è Activation overpotential: 
This is given by the most general form of Butler-Volmer 
equation [46]: 


ap aF. Qf,aF. 

İra = 10,4 [exp ( es re) exp ( a ma) | (50) 
ab cF. F 

ire = ioc fex( b,c pee) exp( fc ree) (51) 


where ap and a are the backward and forward transfer 
coefficients, and ioa and io, are the exchange current densities 
at each electrode, for the calculation of which the following 
experimental correlations are widely used [47,48]: 


: iE 

loa = Ya XH, xtro exp( - ate (52) 
f E 

loc = Yc XO, exp ( = ae) (53) 


with Eacta and Eact, being the anodic and cathodic activation 
energies; a, b and c expressing the concentration dependency; 
and ya and y, being the anodic and cathodic pre-exponential 
coefficients respectively. 


2.5.2. SOEC mode 
If the cell operates as an SOEC, the electrochemical reaction is 
reversed: 


H20 => H2 +} Oo 

The cell voltage is given by Eq. (27), where the open circuit voltage 
of the reverse reaction is to be considered: Vocsogc = — Vocsorc: 
Hence, 


V = Vocsoec Nconc,a — "conc,c ~ Nact,a ~ “Nact,c — ohm <0 (54) 


The overpotentials are positive magnitudes, and their definitions 
slightly differ from those in the SOFC mode, as indicated below: 


© Concentration overpotential: 


_ RT, Px,0.pPHy.ra 
Nconc,a = zp” {Pesci (55) 


RT. V Por 
"conc,c = zp” a (56) 
V Poz,b 


© Activation overpotential: nacta and Nace are given by Eqs. (50) 
and (51) respectively, where ioa and ioc are given by Eqs. (52) 
and (53), since the exchange current density is by definition 
the same for the direct (SOFC) and reverse (SOEC) reactions; 
however the charge transfer coefficients change according to 


Op,a,SOEC = &f,a,SOFC (57) 
@db,c,SOEC = @f,c,SOFC (58) 
@f,a,SOEC = @b,a,SOFC (59) 
Qf, c, SOEC = @b,c,SOFC (60) 


© Ohmic overpotential: nonm is defined in Eq. (43) aS nohm = 
Pera — Perc- IN an SOEC, the current flows from the fuel electrode 


to the air electrode and thus 


Nohm > 0> Pera > Perc = pe = he (61) 


In summary, the electrochemical model consists of a set of 
algebraic equations, namely Eqs. (27), (28), (29), (49), (50) and (51) 
for SOFC mode and Eqs. (27), (55), (56), (49), (50) and (51) for 
SOEC model. However, the transfer coefficients (ap,a;) and the 
exchange current densities (io.a, Îoc) in Eqs. (50) and (51) are cell 
dependent and not measured experimentally. Hence, for comple- 
tely specifying a model, these values are often guessed in a data 
fitting process; in this work, this is presented in Section 5. 


3. Numerical details 


The numerical solution of the mathematical model introduced 
in Section 2 has been addressed by developing an in-house 
algorithm solvable in OpenFOAM-1.5-dev [30] using a finite- 
volume method. In this work we study the steady-state SORFC 
operation, and therefore the time derivatives in the transport 
equations are null. 

To facilitate the implementation of the mathematical model, 
the spatial domain is split into subdomains for each of the five cell 
components; each subdomain is meshed separately resulting in 
five conformal adjacent meshes (Fig. 1), namely: (i) fuelChan- 
nelMesh for the fuel channel domain; (ii) anodeMesh for the fuel 
electrode; (iii) electrolyteMesh is the mesh for the electrolyte 
subdomain; (iv) cathodeMesh for the air electrode; and (v) 
airChannelMesh for the air channel. These numerical grids are 
adjacent to each other; between two adjacent meshes there is only 
one physical boundary but this is represented twice in the 
numerical model, one for each of the two adjacent meshes (these 
are the so-called coupled boundaries). 

The overall algorithm linking the several solution procedures is 
embodied in a main code. Its salient algorithmic steps are shown 
in Fig. 2. First, the computational meshes are created and checked, 
returning a fatal error if they are not face-matching. Then, the 
fields required in each spatial domain are created and their values 
initialised. Next, the radiation model is solved providing the 
radiative heat fluxes to be introduced as a source term in the 
energy-conservation equation for the air electrode and the elec- 
trolyte [16]. Hence, from Eq. (24) the energy-conservation equa- 
tion in the electrolyte, considering surface-to-surface radiation 
and steady state, is 


2 
1 > 
-V-GVI=Z-V- q rad (62) 


where Q aa is non-zero only at the radiative boundaries of the 
electrolyte. Similarly, from Eq. (21) the energy-conservation equa- 
tion in the air electrode, considering surface-to-surface radiation 
and steady-state conditions, is 


a 
v. [TE Cow N o| — V . (af VT) 
Va 


> => > 
= TV * [E CWN o] E dlhaV kj (W N a)]— V. q rad (63) 
where @ „aa is non-zero only at the radiative boundaries of the air 
electrode. 

The electrode model is solved at both electrodes and the 
relevant information, such as velocity, pressure and species con- 
centrations at the coupled boundaries, is transferred to the 
adjacent channel sub-domains. Chemical reaction is modelled 
as a surface reaction; therefore, in Eq. (11) Sx, =0 and the mass 
source/sink is imposed at the reaction wall as a fixed gradient 
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Air channel 


Air electrode 


Electrolyte 


Fuel electrode ERR 


Fuel channel 


Fig. 1. Schematic of the numerical multi-domain. 


SORFC-solver 


{ 


for(startTime, startTime < endTime, startTime +=deltaTime) 
{ 
create: fuelChannelMesh, anodeMesh, electrolyteMesh, cathodeMesh, 
airChannelMesh; 
checkCoupledPatches; 
create: fuelChannelFields, anodeFields, electrolyteFields, cathodeFields, 
airChannelFields; 
fieldsInitialization; 


radiationSolve; 
electrodesSolve; 


electrodeToChannelCoupling; 


channelsSolve; 
channelToElectrodeCoupling; 


electrochemistrySolve; 
electrolyteSolve; 


checkBalances; 
write; 


exit 


Fig. 2. Structure of the segregated solver for the mathematical model of an SORFC. 
boundary condition given by from Eq. (63): 


jens gates V. [TE Cow N o) =y GF YVT) 


eff 
Dap T 


VPar = -erf 2; 


BAa a z5 ee 
=TV. [E CWN o] = Diha V 7 (We N a)I =y: q rad 7 V 
Va Va 


+ 


a 
N ar 1 Par Bo i (64) 


pI i RT; fir Pr 


Kar Kar CE rev + F irrev) (65) 


Ss or . where 
where N ar is given by Faraday's law. Similarly, in Eq. (21) #=0 


and the reaction heat is introduced as boundary heat fluxes. Thus, Grev.c = 9 (66) 


710 

dreva = [AH rr = 0K) + nFVoc]| N uetra -T ral (67) 
irrev,c = — İrc(Nact,c + Ncon,c) (68) 
Girrev,a = — İra(act,a + Ncon,a) (69) 


qrev being non-zero only at the fuel electrode reaction surface, it is 
exothermic if SOFC and endothermic if SOEC; and qirrey being non- 
zero at both electrode reaction walls, it is exothermic at both 
electrodes and for both operation modes (SORFC). These fluxes are 
defined normal to the reaction surfaces. 

Next, the channel model is solved for both the fuel and air 
channels using the SIMPLE velocity—pressure-coupling algorithm 
to solve the mass- and momentum-conservation equations. The 
energy-conservation equation is, at this point, conveniently solved 
for the entire SORFC domain using a coupled solver. Information is 
then exchanged at the channel-electrode coupled boundaries. 

Finally, the electrochemical model equations are solved to calcu- 
late the concentration overpotentials explicitly, and the activation 
overpotentials implicitly using a Newton-Raphson method. Then, the 
value of the electric potential in the ion-conducting phase (electro- 
lyte) at the reaction boundaries is calculated as follows: 


1/2 
SOFC: g= as In Ga Nacte (70) 
rej 

RT Phra 

SOFC: ža = -E In| 5 | +itactat+V (71) 
era JF (2 =) lacta 

RT p 1/2 

SOEC : ere = -JF In (= ) + Nact,c (72) 
2 ,FC 

RT Pu, 0,ra 
SOEC : x a= E In |— V 173 
Pera 2F | Pira | act,a (73) 


These values are set as at the corresponding boundaries of the 
effective electric-potential field in the electrolyte; they are needed for 
the solution of the charge-transfer conservation equation in the 
electrolyte model, which is solved next, thus providing the current 
density distribution in the cell at the preset operating voltage. Mass 
and heat balances are finally checked to ensure conservation at 
convergence. 


4. Experimental data 


The applicability and validity of the model and numerical tool 
is shown in the following sections by comparing the numerical 
results with the data measured by Laguna-Bercero et al. [10] for 
their single, self-supported, hydrogen-fed, microtubular SORFC. 

The experimental test rig is sketched in Fig. 3. It consists of a 
cylindrical furnace, at the centre of which a single microtubular 
SOFC is placed. The furnace has a metal casing and a ceramic 
cylindrical wall with heating resistances. The space between both 
surfaces is filled with a thermally insulating material. The micro- 
tubular SOFC is longer than the furnace, so that both ends are kept 
away from the heated area and economical and simple sealing 
mechanisms can be used. The air electrode, which is shorter than 
the electrolyte/fuel-electrode tube, is located in the middle of the 
heated perimeter to keep the temperature of the cell active area 
close as possible to the furnace preset one. The relevant features of 
the experimental test-rig and its operating conditions are sum- 
marised in Table 4. 

The experimental i-V curves reported by Laguna-Bercero et al. 
[10] are plotted in Fig. 4. We will use in this paper their 


Fuel electrode + Electrolyte 
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Metalic casing 


Heated wall 


Air electrode 


Fig. 3. Schematic of the test rig. 


Table 4 


Thermal insulation 


Reference microtubular SORFC: geometry and operating conditions. 


Microtubular geometry SORFC 

Fuel channel internal radius [10] 1.2mm 

Fuel channel length [10] 10cm 

Anode thickness [10] 400 pm 

Anode length [10] 10cm 

Anode porosity [54] Ea 41.4% 

Anode tortuosity factor [64] Ta 3 

Anode mean pore diameter [54] dpa 1.2 pm 

Anode permeability [64] Boa 7.65e—15 m? 

Electrolyte thickness [10] 20 pm 

Electrolyte length [10] 10cm 

Cathode thickness | 10] 80 pm 

Cathode length [10] 1cm 

Cathode porosity [54] êz 37.1% 

Cathode tortuosity factor [64] Te 5 

Cathode mean pore diameter [54] dpe 0.12 pm 

Cathode permeability” Bo,c 5.5e—17 m? 

Air channel external radius* 1cm 

Air channel length? 7.5 cm 

Solid materials properties SORFC 

Anode thermal conductivity [65] Wea 3Wm`™! K! 

Electrolyte thermal conductivity [65] 2e 2Wm`™! K! 

Electrolyte emissivity* Eo 0.2 

Ionic conductivity constant [65] Ae 85 000 S m7! 

Ionic conductivity constant [65] Be 11 000 K 

Cathode thermal conductivity [65] Ase 4Wm-'!K"! 

Cathode emissivity [66] E 0.4 

Furnace emissivity“ Soven 0.8 

Reference operating conditions SOFC SOEC 

Furnace set temperature [10] Ton 1168 K 1168 K 

Operating pressure [10] Pout 100 000 Pa 100 000 Pa 

Fuel temperature at inlet* Tinet 473K 473 K 

Fuel composition, Hz : H20 : N2[10]°  Xainfuet 15:70:15 15:70:15 
80:20:0 16:20:64 

Fuel flow-rate (at 70%H20, Troom) [10] Qinfuet 0.21min™! 0.21 min~! 

Fuel flow-rate (at 20%H20, T;oom)* Qinjuet 0.1251 min~! 0.1251 min~! 

Air composition, O2 : N2 Kain 0.21:0.79 0.21:0.79 

Air flow-rate (at Troom)* Qinar = 0.11min~' = 0.11 min~! 


* Personal communication. 


P Estimated with the Poiseuille-flow equation [67]. 


© Typical value assumed. 


experimental results for the SOEC mode up to 1.2V since the 
authors subsequently detected electrolyte degradation for higher 
voltages [49]. Also, the experimental open circuit voltage for the 
SOEC at 20% H20(~ 0.87 V) is lower that what it is expected 
according to the Nernst equation ( ~ 0.89 V). Laguna-Bercero et al. 
confirm’ that this mismatch is attributable to an inaccurate 


1 Personal communication. 
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SOFC-SOEC, T = 1168 K 


SOEC,,,: PH20 = 0.7 atm 
SOFC,,5: PH20 = 0.7 atm 
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° 
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SOECexp: PH20 = 0.2 atm 
SOFCexp: PH20 = 0.2 atm 


DO 


Fig. 4. Experimental characteristic curves of the microtubular SORFC of Laguna-Bercero et al. [10]. 


performance of the humidifier. The actual composition of the 
feeding gases has been instead estimated to be 14% H2, 30% H2O 
and 56% N3. This is therefore the composition used to calculate the 
numerical curve labelled “SOEC,im : pH20 = 0.2 atm”. 


5. Fitting process 


In this section we describe the fitting process employed to 
assign values to the electrochemical parameters of a particular cell 
prior to using the model for its simulation. 

The need of this fitting process stems from the multiscale 
nature of the several phenomena in the SORFC. Thus CFD methods 
usually resolve the macroscale phenomena, while the microscale 
processes are modelled and included as source terms in the 
transport equations. This is the case of the electrochemical reac- 
tion rate, given by the kinetic equations of Butler-Volmer, 
Eqs. (50) and (51). These involve the following electrochemical 
Parameters: dba, Ata, Abc, Ate, loa, toc, the values of which are 
difficult to acquire experimentally. 

In the CFD modelling literature, the charge transfer coefficients 
are often set to 0.5 both for SOFC [15,50,12] and SOEC [17,19,51]. 
This corresponds to a simple, single-electron-exchange reaction- 
mechanism [46], despite the fact that hydrogen-oxidation reaction 
involves 2 electrons. Nevertheless, this simplification is generally 
accepted because of its ease of use and apparently good results. 
The present authors did successfully use this assumption in a 
previous modelling study for SOFCs [16]. However, we find that 
this value may be inappropriate when modelling the SOFC and 
SOEC modes simultaneously, as shown below. 

The estimation of the exchange current density is typically the 
other key issue in the parameter fitting process. In the simplest 
approach among those used in the literature, the exchange current 
density is assumed constant and its value is fitted, as in [20]. In more 
comprehensive models, experimental correlations are used to account 
for the exchange-current-density dependence on temperature and 
species concentrations. Eqs. (52) and (53) represent typical correla- 
tions for the exchange current densities of both electrodes, where ya, Ye 


are the fitted parameters and the values a, b, c significantly vary from 
one correlation to another due to their strong dependence on the cell 
materials and microstructure. None of the CFD models reviewed by 
the authors using one of these empirical correlations from the 
literature justify the accuracy of the chosen correlation. 

In the following subsections, we present the parameter fitting 
process performed in this work to estimate the values of the 
charge-transfer coefficients and the exchange current densities. 
The aim is to provide more reliable fitted values for SOFC and SOEC 
that can be used in a quantitative modelling tool. 

5.1. Estimation of the charge transfer coefficients 

The charge transfer coefficients (apa, fa, @b,c, fc) are usually set 
to 0.5 (i.e. single step, single electron transfer) [15,50,12,17,19,51], 
for the reasons, and despite the shortcomings, referred to above; 
therefore, for most models the fitted parameters are just the 
exchange current densities. 

Because of our previous successful experience with these 
values for SOFCs [16], we have in the first instance attempted to 
use them for the present SOFC-SOEC model. Fig. 5 shows the 
experimental data of Laguna-Bercero et al. [10] together with the 
numerical results obtained when taking ap a = ata = Apc = afc = 0.5 
(see “case-0” in Table 5); the numerical value of i is calculated over 
the air—-electrode active area as 


=E f Tre . Tre dArc 


From the figure, it is apparent that taking all the charge transfer 
coefficients as 0.5 seems to work for the SOFC mode (as in [16]) 
but fails to model the SOEC operation, as indicated by the arrows 
in magenta. According to the Butler-Volmer equation (Eq. (50)), 
the magnitude of the current density depends on the forward and 
backward charge transfer coefficients (apa, afa, @b,c, Mc) regardless 
of whether the operating mode is SOFC or SOEC (forward or 
backward). This is not satisfied by the numerical results shown in 
Fig. 5. Therefore, in the present study, we reject the hypothesis 
Qba = Ofa = Ap, = @c =0.5 and we estimate the value of each 
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SOFC-SOEC, T = 1168 K: Oy g = Ofa = Uy, = Ofo = 0.5 
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Fig. 5. Characteristic curves of a microtubular SORFC. Experimental data (open symbols) for the cell in Table 4 [10] and simulation results (solid symbols) for different values 
Of aba = fa = Apc = arc (see Table 5). Arrows highlight discrepancies between measurements and modelling. (For interpretation of the references to color in this figure 


caption, the reader is referred to the web version of this paper.) 


Table 5 


Electrochemical parameters for both operating modes and for the different cases studied. 


Parameter SOFC SOEC SOFC SOEC SOFC 
Ref. Ref. Case-0 Case-0 Case-1 
aba 1 0.1 0.5 0.5 1 
aba 0.1 1 0.5 0.5 0.1 
Eicta 120 120 120 120 120 
Ya 2.5e9 2.5e9 2.5e9 2.5e9 2.5e9 
a 0.11 0.11 0.11 0.11 1 
b 0.67 0.67 0.67 0.67 —0.5 
fa 0.25 0.25 0.25 0.25 0.25 
abc 0.4 0.1 0.5 0.5 0.4 
aic 0.1 0.4 0.5 0.5 0.1 
E;cte 120 120 120 120 120 
Ye 2.6e9 2.6e9 2.6e9 2.6e9 2.6e9 


charge transfer coefficient that best fits the experimental data. The 
resulting values are reported in Table 6; for SOFC mode these 
are: aba = 1, ata = 0.1, abc = 0.4, atc =0.1. For SOEC mode the 
forward and backward directions are inverted and hence: ap a= 
0.1, ata = 1, Apc = 0.1, atc = 0.4. 

These results may be quantitatively interpreted according to 
the theory of the Butler-Volmer equation for multistep reactions 
[46], which states that 


Qb + OF = us (75) 
V 

with n being the number of electrons exchanged in the electro- 
chemical reaction and v the number of occurrences of the rate- 
determining step in the overall reaction. For a hydrogen fed SORFC, 
such the one studied in this work, n=2. Hence, from Eq. (75) ve = 4 
and va = 1.8, meaning that the limiting step in the air electrode 
occurs four times and the one in the fuel electrode happens almost 
twice per each overall reaction. 

These values agree qualitatively with those reported by 
Zhu and Kee [45] (SOFC mode: apq=1.5, ata = 0.5, abc = 0.75, 
ac = 0.25) for a button SOFC with similar materials as those used 


SOEC 
Case- 


SOFC SOEC SOFC SOEC SOFC SOEC 
1 Case-2 Case-2 Case-3 Case-3 Case-4 Case-4 
1 0.1 1 0.1 1 0.1 
0.1 1 0.1 1 0.1 1 
120 120 120 120 120 120 
1e10 1e10 2.5e9 2.5e9 5e9 5e9 
1 1 1 1 1 1 
-0.5 —0.5 —0.5 —0.5 —0.5 —0.5 
025 0.25 0.25 0.25 0.25 0.25 
0.4 0.1 0.4 0.1 0.4 0.1 
0.1 0.4 0.1 0.4 0.1 0.4 
120 120 120 120 120 120 
2.6e9 2.6e9 1e10 1e10 1.8e9 1.8e9 
Table 6 
Reference microtubular SORFC: electrochemical parameters for both operating 
modes. 
Electrochemical parameters SOFC SOEC 
Fuel transfer coefficients* Aba, Oa 1,01 0.1, 1 
Fuel activation energy [48] Eata 120kJmol~! 120 kJ mol! 
Fuel pre-exponential coefficient? Ya 2.5e9 A m~? 2.5e9 A m~? 
Exponent [47] a 0.11 0.11 
Exponent [47] b 0.67 0.67 
Exponent [48] c 0.25 0.25 
Air transfer coefficients? Apc, Ae 0.4, 0.1 0.1, 0.4 
Air activation energy [48] Eacte 120 kJ mol™! 120 kJ mol`! 
Air pre-exponential coefficient? Yc 2.6e9 A m~? 2.6e9 A m~? 


* Fitted parameters. 


by Laguna-Bercero et al. [10], i.e. Ni-YSZ, YSZ, LSM. Our results also 
agree with the experimental data of Kawada et al. [23] for a Ni-YSZ 
SOFC anode (SOFC mode: apa = 2, ata = 1) and Kenney and Karan 
[24] for a LSM/YSZ SOFC cathode (ap, = 1.5, a¢¢ = 0.5). The agree- 
ment is particularly revealing in two respects: (i) the forward 
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charge transfer coefficient is smaller than the backward one for A corollary of this research is that the often-employed assump- 

both electrodes; and (ii) the charge transfer coefficients are larger tion apa = Aq = Ab, = afc = 0.5 fails when properly validated with 

in the fuel electrode than in the air one. a complete set of experimental data for the cell operating in direct 
a 
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Fig. 6. Characteristic curves of a microtubular SORFC. Experimental data (open symbols) for the cell in Table 4 [10] and simulation results (solid symbols, see Table 5): 
(a) varying ipa; (b) varying ioc; (c) varying simultaneously ioa and ioc. (For interpretation of the references to color in this figure caption, the reader is referred to the web 
version of this paper.) 


714 


and reverse mode (SORFC). If the model is only validated with 
SOFC or SOEC data, the error made by this assumption will be 
compensated by the data fitting, ie. the exchange current densi- 
ties, and this in turn can over/under-estimate the activation 
overpotentials. 


5.2. Estimation of the exchange current density 


The exchange current density is often estimated using empirical 
correlations such Eqs. (52) and (53), where a, b, c are set by the 
correlation and then ya, yc are the only fitted parameters. In the works 
reviewed, c is generally set to 0.25 while a and b vary significantly; for 
instance Yamamura's model specifies a=1 and b=—0O.5 (used in 
{48,28,52]), Mogensen's model (used in [48,51,17]) chooses a=b=1 
and recent experimental works report lower values for a and inter- 
mediate values for b (a=0.11, b=0.67) [47,53]. 

Yamamura's values were used in our previous work [16]. 
Validation was restricted to a single set of SOFC data, and an 
acceptable agreement was found. Thus these values were used in a 
first attempt in the present work to simulate a similar reversible 
cell for two different fuel compositions. Given a=1, b=—0.5 
(Yamamura) and c=0.25, the value of ya and yc is guessed in a 
trial and error method aiming at the best possible agreement 
between numerical and experimental i-V curves. Representative 
results from this fitting process are plotted in Fig. 6. 

Fig. 6(a) shows the experimental data [10] (open symbols) and 
two iterations of the fitting process, labelled “case-1” and “case-2” 
and specified in Table 5. As indicated in the figure by magenta 
arrows, “case-1” (solid line with symbols) underpredicts the mean 
current density of the SORFC with 70% H20 and the SOEC with 20% 
H20; on the contrary, it overpredicts the performance when 
running the cell as an SOFC with 20% H20. In order to get a better 
fitting, another value of ya is assumed in “case-2” (dashed line with 
symbols). The results of “case-2” improve those of “case-1” 
reducing the performance underestimation of the latter and 
providing a good fit for the whole SORFC operating range under 
conditions of 70% H20. However, as it is pointed out with dashed 
arrows in cyan, “case-2” overpredicts the current density for both 
SOFC and SOEC operation with 20% H20. 
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Fig. 6(a) shows that when ya is modified to improve the fit for 
curves corresponding to high steam concentrations in the fuel, 
then the fit of the curves at low steam concentrations worsens. 

Similar to Fig. 6(a), Fig. 6(b) shows the results of the fitting process 
as the value of y, is changed; this is increased from “case-1” (solid line 
with symbols) to “case-3” (dashed line with symbols). It is found that 
“case-3”(see cyan dashed arrows): (i) provides a good fit for the 70%- 
H20 SOFC and the 20%-H20 SOEC curves; (ii) slightly improves the 
agreement between the 70%-H20 SOEC curves; and (iii) worsens the 
agreement of the 20%-H20 SOFC curves. 

Fig. 6(b) shows that varying the value of yc is not enough to get 
simultaneously a good fit between numerical and experimental 
data for the four curves considered. 

Finally, Fig. 6(c) shows the results when ya and y, are modified 
simultaneously. Here “case-4” (dashed line with symbols) is obtained 
with a lower value of ye and a higher value of ya than those used in 
“case-1” (solid line with symbols). The 20%-H20 SORFC performance is 
well predicted in “case-4” but it fails in the prediction of the 70 %-H20 
SORFC, as indicated by the dashed cyan arrows. 

After several iterations like those represented by “case-1”, “case-2”, 
“case-3” and “case-4”, we conclude that there is no pair of values for ya 
and y, that makes the numerical and experimental data agree simu- 
Itaneously in all four situations studied. However, some combinations 
provide a good fit of the SOFC and SOEC curves under the same 
operating conditions, e.g. “case-2” reproduces the performance of the 
70 %-H2O SORFC and “case-4” that of the 20%-H,0 SORFC. Therefore, 
we conclude that the model can describe the SORFC performance but 
not its dependence on the fuel concentration. This conclusion points to 
Yamamura's values as the cause of this underperformance. 

The work of Bieberle et al. [53], recently supported by Utz et al. 
[47], was then used to replace Yamamura's correlation. A fitting 
process similar to that presented above was conducted using 
Bieberle's values of a=0.11 and b=0.67, resulting in a pair of 
values for y, and y, that provide good fit of the experimental data 
both for low and high steam concentrations and in both SOFC and 
SOEC operations. The best fit is that shown in Fig. 7 for the 
electrochemical values reported in Table 6. 

To summarise this section, we have shown that the results of a 
macroscale CFD tool strongly depend on the correlation chosen to 
model the dependence of the exchange current density with the gas 
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Fig. 7. Characteristic curves of a microtubular SORFC. Experimental data (open symbols) [10] and simulation results (solid symbols) for the cell in Tables 4 and 6. 
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Fig. 8. Mean overpotentials for all the simulations plotted in Fig. 7. 
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Fig. 9. Mean exchange current density over the electrodes active area for both electrodes (ioa, io.) and for all the simulations plotted in Fig. 7, together with the cell average 


current exchange density (io cet) = (io. +toc)/2) and a reference experimental value [55]. 


composition at the microscale. The authors recommend the use of 
(at least) a couple of experimental data sets with different reactant 
compositions to check whether a correlation from the literature, 
expressing this dependence, truly applies to a cell being modelled. 


6. Validation 


The validation of the mathematical model and numerical tool 
presented in this work is shown in Fig. 7, where the experimental 


data for a single, self-supported, hydrogen-fed, microtubular solid 
oxide regenerative fuel cell [10] is compared with the numerical 
results obtained after the fitting process presented in the former 
section. The values of the fitted parameters are reported in Table 6, 
namely: (i) the pre-exponential coefficients (ya, yc), the exponents 
(a, b) in Eq. (52), and the backward and forward transfer coeffi- 
cients (dba Qf a, b,c; afc). 

As shown in Fig. 7, numerical and experimental curves are in 
good agreement. However, the applicability of such numerical 
tool must be further verified to check the reliability of the values 
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estimated for the parameters. In the literature, however, this 
additional validation is rarely performed, probably because of 
the lack of complementary experimental data to validate against 
(e.g. Electrochemical Impedance Spectroscopy, EIS, analysis). 

Fig. 8 shows the distribution of the voltage losses obtained 
numerically for the given values of the fitted parameters. It is 
found that the main voltage drop is due to the activation losses, as 
it may be presumed from experimental data [54,49]. Moreover, 
Fig. 9 shows the distribution of the calculated mean exchange 
current densities, which are in the same order of magnitude as the 
value estimated from experimental EIS data (ig ce + 0.5 A cm~?) 
for a similar cell running as SOFC under slightly different operating 
conditions (V=0.95 V, T=1168 K, H2: H20=4.85:3) [55]. It is thus 
shown that the values of the fitted parameters not only provide a 
good accord between the numerical and experimental character- 
istic curves, but also reliable results of the internal physics of 
the cell. 


7. Conclusions 


In this work, we have presented a comprehensive numerical 
tool for the simulation of a solid oxide regenerative fuel cell. The 
tool models the following SORFC phenomena: (a) mass transport 
through channels and porous solids; (b) momentum transport 
through channels and porous solids; (c) species transfer through 
channels and porous media (convection, ordinary diffusion, Knud- 
sen diffusion); (d) heat transfer in gases and solids by conduction, 
convection and surface-to-surface radiation; (e) charge transport 
through an impervious solid; and (f) reversible electrochemistry. 

The validation of such numerical tool involves the fitting of 
the unknown parameters in the electrochemical model, namely: 
(i) the pre-exponential coefficients (ya, yc), the exponents (a, b) in 
Eq. (52), and the backward and forward transfer coefficients 
(aba, Qf,a, @b,c, Oc) 

A literature review regarding this fitting process is also pre- 
sented. It leads to the conclusion that this process often over- 
looked in the literature, where it is (if at all mentioned) not clearly 
explained or sufficiently justified. Thus, the fitting is here 
described in detail. Two corollaries emerges in this research. 

First, the often-employed assumption ap a= afa = abc = ac = 0.5 
fails when properly validated with a complete set of experimental data 
for the cell operating in direct and reverse mode (SORFC). Our results 
suggest that the forward charge transfer coefficient is smaller than the 
backward one for both electrodes; and that the charge transfer 
coefficients are larger in the fuel electrode than in the air one. 

Second, the relevance of exchange-current-density modelling is 
also highlighted in this study; we have shown that an accurate 
modelling of the dependence of the exchange current density on 
the reactant concentration is paramount for a reliable prediction of 
the overall cell performance. 

The validity of the numerical results provided after the fitting 
process for some of the uncertain parameters has been shown by 
comparison with experimental data for an anode-supported 
microtubular SORFC. Good agreement is found not only between 
the numerical and experimental i-V curves, but also between 
numerical and experimental electrochemical parameters, such as 
activation overpotentials and exchange current densities, which 
supports the accuracy of the fitted values. 
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Appendix A. Thermophysical properties 
A.1. Viscosity, u 


The semi-empirical formula of Wilke is used to estimate the 
viscosity of a multicomponent fluid [68]: 


Xah 
= a A.1 
i a il ii 
where 
1/2 1/442 
y= [1+ (a/Ug) (W/W) 1] (A.2) 


(8+8W,/W,)'/? 


The Sutherland model accounts for the temperature depen- 
dence of the viscosity , of each species a: 


I Topa oF Siia T 3/2 
Fa = Mog T + Sja Towa 


(A.3) 


where the Sutherland-law parameters (79, To,,S,,) are tabulated for 
the most common gases in [68]. 


A.2. Binary diffusion coefficient, Dap 


The binary-diffusion coefficients are modelled according to the 
empirical correlation given by Wilke [69]: 
3Wa +Ws 
ZW, W; 
O2,2dP 


Dap = 2.628 x 107? (AA) 


where the binary diffusion coefficient Dag is given in square 
centimetres per second; the temperature, T, in Kelvin; the mole- 
cular weight of species a, Wa, in kilogram per kilomole; the 
pressure, p, in atmospheres. The average collision diameter, in 
angstroms, is calculated as 


Oat Op 
2 


where the value of o, (or og) is tabulated for the most common 
species in [70]. Qp is the dimensionless collision integral in the 
Lennard-Jones 12-6 potential model: 


Cup = (A.5) 


1.06036 0.19300 | 1.03587 _ 1.76474 AG 
D= 770.15610 " 9(0.47635T) * ẹ(1.52996T) ` e(3.89411T) (A.6) 
Here, 
r- (A7) 
/EaEp 


where kg is the Boltzmann constant and e, is the characteristic 
Lennard-Jones energy of species a. For the most common gases the 
value of ¢, is tabulated in [70]. 

To convert the units of the binary diffusion coefficient and the 
pressure into the SI system (m° s71, Pa), Eq. (A.4) is modified as 
follows: 


2.628 x 1073, /T3Wat Wo 
Dap = 10.1325 V we 


2 
Tgp DP 


(A.8) 


where Gap is the only parameter that is not expressed in SI system 
units, as it remains in angstroms. 
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The diffusion coefficient of a species a in a multicomponent 
mixture can be calculated as 


=x; 


Sonal ce 
BAa Dap 


Xa being the molar fraction of species a and D,, being the binary 
diffusion coefficient defined above. 


Dam (A.9) 


A.3. Knudsen diffusion coefficient, Dga 


The Kinetic Theory of Gases defines the Knudsen diffusion 
coefficient as [71,72] 


dp |8RT 
Ka "3 aW, 


(A.10) 


where Dx, is given in SI units if the temperature is expressed in 
Kelvin, the molecular weight in kilograms per kilomole, the mean 
pore diameter (dp) in meters and the ideal gas law constant in 
Joule per kilomole per Kelvin. 


AA. Thermal conductivity, A 


Similar to viscosity (Section A.1), the thermal conductivity of 
the fluid is given by Wilke's formula [68]: 


Xad 
= | A.11 
4 D T] a 
where 
1/2 1/442 
bap = LTO) W/W) "1 (A.12) 


(8+8W,/W,)'/ 


The thermal conductivity of each species a, Aq, is then given by the 
Sutherland model: 


Tora t+Sia ( T ) 3/2 


4a = doa T+Sia \Toia 


(A.13) 


where the Sutherland-law parameters (20, To, S,) are tabulated for 
the most common gases in [68]. 


A.5. Specific heat at constant pressure, C, 


The specific heat at constant pressure of a multicomponent gas 
is calculated as 


Cp = Z YaCpa (A.14) 


where Cpa, the specific heat at constant pressure of species a, is 

provided by the JANAF thermochemical tables: 
R 2 3 4 

Cpa = (Mat Geel +03aT +44,T +5 1 } (A.15) 


The parameters aia, 2, 43a, A4a, A5a are the JANAF constants [73]. 


A.6. Sensible enthalpy, h 


The sensible enthalpy of a multicomponent fluid is given by 


h=Laha) (A.16) 
Va 
where h, is the sensible enthalpy of each component: 
T 
ha= | Coa AT (A.17) 
Tref 


The reference temperature is for convenience set to zero, 
Tef = 0 K. Hence, from Eqs. (A.17) and (A15): 


ize R (ara T+ Ber? + Ber? 4 Sars + Bers) (A.18) 


WwW, 2 3 4 5 
A.7. Ionic conductivity, 6 
The electrode is the only ion-conductor SOFC-component. 


Theory and experiments suggest for ion-conducting electrolyte- 
materials the following equation for the ionic conductivity: 


Ge =Acexp ( -") 


where the constants Ae and Be are tabulated for most common 
SOFC electrolytes in [65]. 


(A.19) 
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